       !=======================================================================
       !> @file modules.f90
       !> @brief Global functions and subroutines module
       !> @authors Ary Rodriguez-Gonzalez & Alejandro Esquivel
       !> @date 23/Jun/2011
       !
       ! Copyright (c) 2011 Ary Rodriguez-Gonzalez & Alejandro Esquivel
       !
       ! This file is part of AGA-V1.
       !
       ! AGA-V1 is free software; you can redistribute it and/or modify
       ! it under the terms of the GNU General Public License as published by
       ! the Free Software Foundation; either version 3 of the License, or
       ! (at your option) any later version.
       !
       ! This program is distributed in the hope that it will be useful,
       ! but WITHOUT ANY WARRANTY; without even the implied warranty of
       ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
       ! GNU General Public License for more details.
       ! You should have received a copy of the GNU General Public License
       ! along with this program.  If not, see http://www.gnu.org/licenses/.
       !=======================================================================
       !> @brief module mpi_module
       !> @details This module contains declarations needed to implement the
       !! Message Passing Interface (MPI) paralelization.
       !-----------------------------------------------------------------------
       module mpi_module

               implicit none

#ifdef MPIP

               include "mpif.h"
       
       !> @brief rank of the master node (who is in charge of the I/O)
       integer, parameter :: master=0
       
       !> @brief rank of each processor (if MPI is disabled has a value of 0)
       integer :: rank
       !> @brief total number of processors (if MPI is disabled has a value of 1)
       integer :: np
       integer :: err                                          !< @brief MPI error variable
       integer :: status(MPI_STATUS_SIZE)      !< @brief MPI status variable
       
#else
       
       integer, parameter :: np=1, master=0, rank=0
       
#endif
       
       end module mpi_module
       !-----------------------------------------------------------------------
       !> @brief module func
       !! @details This module contains global functions and subroutines used
       !! by AGA-V1
       module func
       contains
       !-----------------------------------------------------------------------
       !> @brief Updates the size of the hyper-boxes
       !> @param[in] npar : number of free parameters (parameter space
       !! dimension, denoted by "n" in the paper)
       !> @param[in] niter : number of iteration counter
       !> @param[in] nind : number of individials in the entire population
       !> @param[in] deltaconst : if equal to 1 the reduction of the hyperboxes
       !! is at a constant rate (see reduc), otherwise hyperboxes are reduced 
       !! with the standard deviation of each parameter.
       !> @param[in] reduc : reducing factor, if a constant speed is chosen the
       !! hyper-box is reduced according to: 
       !! @f$\Delta_{i,niter}=\Delta_{i,niter-1}*reduc^{niter}@f$
       !> @param[in] arrmain(npar,nind) : Entire population of "nind" individuals, 
       !! each with "npar" parameters.
       !> @param[in] delta0(npar) : Original size of the hyperboxes, 
       !! array of @f$\Delta_{i,1}@f$.
       !> @param[out] delta(npar) : array of @f$\Delta_{i,niter}@f$ (niter 
       !! is constant, thus is only a vector w/npar entries). 
       !!@f$\Delta_{i,niter}@f$ is the size of the hyper-box corresponding to 
       !! the "i-th" parameter
         subroutine deltabox(npar,niter,nind,deltaconst,reduc,arrmain,delta0,delta)
           implicit none
           real,dimension(npar), intent(in)  :: delta0
           real,dimension(npar), intent(out) :: delta
           real,    intent(in),  dimension(npar,nind)       :: arrmain
           integer, intent(in):: npar,niter,deltaconst,nind
           real, intent(in) :: reduc
           real :: mean
           integer :: i
           
           if (deltaconst .eq. 1 .or. niter .eq. 1) then
              do i=1,npar
                 delta(i) = delta0(i)*reduc**niter
              enddo
           else
              do i=1,npar
                 delta(i)=sigm(arrmain(i,:),nind)
              enddo
           endif
       
         end subroutine deltabox
       !-----------------------------------------------------------------------
       !> @brief Computes the mean value
       !> @param[in] N : number of data points
       !> @param[in] stat(N) : vector from which the mean is to be obtained
       !> @return Arithmetic mean value of stat(N)
          real function med(stat,N)
            implicit none
            integer, intent(in)   :: N
            real, dimension(N), intent(in) :: stat
            integer :: i
            !
            med=0.
            do i=1,N
               med=med+stat(i)
            enddo
            !
            med=med/N
            !
           return
           !
         end function med
       !-----------------------------------------------------------------------
       !> @brief Computes the standard deviation
       !> @param[in] N : number of data points
       !> @param[in] stat(N) : vector from which the standard deviaiton is to 
       !! be obtained
       !> @return standard deviation of stat(N)
          real function sigm(stat,N)
            implicit none
            integer, intent(in)   :: N
            real, intent(in), dimension(N) :: stat
            real :: mean
            integer :: i
            !
            mean=med(stat,N)
            !
            sigm=0.
            do i=1,N
               sigm=sigm+(stat(i)-mean)**2
            enddo
            !
            sigm=sqrt(sigm/(N))
            !
            return
            !
          end function sigm
       !-----------------------------------------------------------------------
       !> @brief Asexual reproduction routine
       !> @param[in] npar : number of free parameters
       !> @param[in] gensize : size of the new generation from each parent.
       !! Since the parent is kept gensize= num of sons + 1
       !> @param[in] nmig : turns on/off migration. If its value is 0, 
       !! migration is disabled, otherwise is enabled
       !> @param[in] ngaussdev : if value equal to 1, offsprinf is obtained
       !! with a Gaussian distribution, otherwise they are obtained with
       !! uniform random deviates.
       !> @param[in] delta(npar) : Size of the parent hyper-box.
       !> @param[in] qpar(npar) : parent parameters (its location in the 
       !! npar-dimensional parameter space
       !> @param[in] qbounds(npar) : bounds of the initial universe
       !> @param[out] qfam(npar,gensize) : A new family (progeny+parent)
         subroutine spawn(npar,gensize,nmig,ngaussdev,delta,qpar,qbounds,qfam)
           implicit none
           integer, intent(in) :: npar,gensize,nmig,ngaussdev
           real, dimension(npar), intent(in):: delta, qpar
           real, dimension(2,npar), intent(in) ::qbounds
           real, dimension(npar,gensize) :: qfam
           integer :: i,j,rand_seed
           real :: ran1,ran
           !
           do i=1,npar
              qfam(i,1)=qpar(i)
           enddo
       
               do j=2,gensize
                       i=1
                       do while (i.le.npar)
                       if (ngaussdev .eq. 1) then
                           ran1=gaussdev(rand_seed)
                       else    
                    call random_number(ran)
                    ran1=2.*ran-1.
                               endif
                               
                               qfam(i,j)=qfam(i,1)+delta(i)*ran1
                 
                               if (nmig.ne.0) then
                                       i=i+1
                               else if ((qfam(i,j) .ge. qbounds(1,i)).and.                 &
                                        (qfam(i,j) .le. qbounds(2,i)) ) then
                                       i=i+1
                               endif
               enddo 
               enddo   
       
       return
       end subroutine spawn
       !-----------------------------------------------------------------------
       !> @brief Evaluates convergence criterion
       !> @param[in] npar : Number of free parameters
       !> @param[in] q(npar) : Individual data
       !> @param[in] delta(npar) : Window size of the individual
       !> @param[in] beta : Tolerance value
       !> @param[out] ncrit : Return value=1 if criterion met, otherwise is not
       !! modified
       subroutine criterion(npar,q,delta,beta,ncrit)
         implicit none
         integer, intent(in)  :: npar
         real, intent(in), dimension(npar) :: q, delta
         real, intent(in) :: beta
         logical, dimension(npar) :: crit
         integer, intent(out) :: ncrit
         integer :: i
       
         crit(:)=.False.
         do i=1,npar
            if (abs(delta(i)/q(i)) .le. beta) crit(i)=.true.
         enddo
       
         if (all(crit(:)).eqv. .true.) ncrit=1
       
         return
       end subroutine criterion
       !-----------------------------------------------------------------------
       !> @brief Generates new realization from experimental data
       !> @param[in] ndata : number of experimental data points
       !> @param[in,out] yobs(ndata) : experimental data
       !> @param[in] ysigma(ndata) : associated uncertainties of yobs
       !> @todo Hay que implementar dispersion Gaussiana y checar por que no se 
       !! guarda yobs0
       subroutine randdata(ndata, yobs, ysigma)
         implicit none
         integer, intent(in) :: ndata
         real, intent(inout), dimension(ndata) :: yobs
         real, intent(in), dimension(ndata) :: ysigma
         integer :: i
         real :: ran1,ran
         !
         do i=1,ndata
            call random_number(ran)
            ran1 = 2.*ran-1.
            yobs(i)= yobs(i)+ysigma(i)*ran1
         end do
       !
       return
       !
       end subroutine randdata
       !-----------------------------------------------------------------------
       !> @brief Generates pseudo-random (Gaussian distributtion)
       !> @details Modified from "Numerical Recipes in Fortran 77, second ed."
       !! (Press et al., Cambridge University Press)
       !> @param[in] idum : dummy variable (required)
       !> @return Pseudo-random number, drawn from a normal distribution
       !! (Gaussian with mean=0 and standard deviation=1)
       function gaussdev(idum)
         implicit none
         real :: gaussdev
         !  
         integer :: iset,idum
         real :: fac,gset,rsq,v1,v2,ran
         save iset,gset
         data iset/0/
         !
         rsq=0.
         if (idum .lt. 0.) iset=0
         if (iset .eq. 0) then
            do while (rsq .gt. 1.0 .or. rsq .eq. 0.)
               call random_number(ran)
               v1=2.*ran-1.
               call random_number(ran)
               v2=2.*ran-1.
               rsq=v1**2+v2**2
            enddo
            fac=sqrt(-2.*log10(rsq)/rsq)
            gset=v1*fac
            gaussdev=v2*fac
            iset=1
         else
            gaussdev=gset
            iset=0.
         endif
         !
         return
       end function gaussdev
       !-----------------------------------------------------------------------
       !> @brief Select the fittest individuals (nparents number of future
       !! parents) from the entire population. Based in the smallest value of
       !! the merit function (smaller number=fittest, e.g. @f$\chi^2@f$)
       !> @param[in] npar : numbr of free parameters
       !> @param[in] nind : number of individuals in the entire population
       !> @param[in] nparents: number of parents, total of fittest individuals
       !! to survive for the next generation
       !> @param[in] arrmain(npar,nind) : entire population
       !> @param[in,out] fitness(nind) : fitness level of each individual, 
       !! given by the merit function provided b the user in "tarfunc.f90".
       !> @param[out] arrpars(npar,nparents) : array of the fittest individuals,
       !! which will be parents in the following generation.
       subroutine pickbest(npar,nind,nparents,arrmain,fitness,arrpars)
         implicit none
         integer, intent(in) :: npar,nparents,nind
         real,    intent(in),  dimension(npar,nind)       :: arrmain
         real,    intent(inout),  dimension(nind)         :: fitness
         real                  ,  dimension(nind)         :: fit2
         real,    intent(out), dimension(npar,nparents):: arrpars
         integer, dimension(nind) :: index
         integer :: i
         
         call indexx(nind,fitness,index)
         do i=1,nparents
            arrpars(:,i)=arrmain(:,index(i))
            fit2(i)=fitness(index(i))
         end do
         
         fitness(1:nparents)=fit2(1:nparents)
         
       contains
       !-----------------------------------------------------------------------
       !> @brief indexing routine
       !> @details Modified from "Numerical Recipes in Fortran 77, second ed."
       !! (Press et al., Cambridge University Press)
       !> @param[in] n: number of elements in the list
       !> @param[in] arr(n) : array to be indexed
       !> @param[in] indx(n) : index list, such that arr(indx(j)) is in
       !> ascending order for j=1, 2, 3 ...n.  arr and n are not changed
         subroutine indexx(n,arr,indx)
           implicit none
           integer, intent(in)   :: n
           real, intent(in),  dimension(n) :: arr
           integer, intent(out), dimension(n) :: indx
           integer :: i,indxt,ir,itemp,j,jstack,k,l
           real    :: a
           integer, parameter :: m=15, nstack=50
           integer, dimension(nstack) :: istack
           !
           do j=1,n
              indx(j)=j
           end do
           jstack=0
           l=1
           ir=n
           do
              if(ir-l.lt.m)then
                 do j=l+1,ir
                    indxt=indx(j)
                    a=arr(indxt)
                    do i=j-1,l,-1
                       if(arr(indx(i)).le.a) exit
                       indx(i+1)=indx(i)
                    end do
                    indx(i+1)=indxt
                 end do
                 if(jstack.eq.0) RETURN
                 ir=istack(jstack)
                 l=istack(jstack-1)
                 jstack=jstack-2
              else
                 k=(l+ir)/2
                 itemp=indx(k)
                 indx(k)=indx(l+1)
                 indx(l+1)=itemp
                 if (arr(indx(l)).gt.arr(indx(ir))) then
                    itemp=indx(l)
                    indx(l)=indx(ir)
                    indx(ir)=itemp
                 end if
                 if (arr(indx(l+1)).gt.arr(indx(ir))) then
                    itemp=indx(l+1)
                    indx(l+1)=indx(ir)
                    indx(ir)=itemp
                 end if
                 if (arr(indx(l)).gt.arr(indx(l+1))) then
                    itemp=indx(l)
                    indx(l)=indx(l+1)
                    indx(l+1)=itemp
                 end if
                 i=l+1
                 j=ir
                 indxt=indx(l+1)
                 a=arr(indxt)
                 do
                    do
                       i=i+1
                       if(arr(indx(i)).ge.a) exit
                    end do
                    do
                       j=j-1
                       if(arr(indx(j)).le.a) exit
                    end do
                    if(j.lt.i) exit   
                    itemp=indx(i)
                    indx(i)=indx(j)
                    indx(j)=itemp
                 end do
                 indx(l+1)=indx(j)
                 indx(j)=indxt
                 jstack=jstack+2
                 if(jstack.gt.nstack) pause 'NSTACK to small to indexx'
                 if(ir-i+1.ge.j-l)then
                    istack(jstack)=ir
                    istack(jstack-1)=i
                    ir=j-1
                 else
                    istack(jstack)=j-1
                    istack(jstack-1)=l
                    l=i
                 end if
              end if
           end do
           !
         end subroutine indexx
       
       end subroutine pickbest
       !-----------------------------------------------------------------------
       !> @todo Preguntarle a ary que chingaos hace esta rutina
       subroutine gaussian(ngaus,i,arrmain)
       implicit none
       integer :: ng,ng2
       real, dimension(3) :: aux
       integer, intent(in) :: i,ngaus
       real, intent(inout), dimension(:,:) :: arrmain
       !
       do ng=ngaus-1,1,-1
          do ng2=0,ng-1
             !
             if(arrmain(3*ng2+2,i) .gt. arrmain(3*ng2+5,i)) then
                aux(1:3)=arrmain(3*ng2+1:3*ng2+3,i)
                arrmain(3*ng2+1:3*ng2+3,i)=arrmain(3*ng2+4:3*ng2+6,i)
                arrmain(3*ng2+4:3*ng2+6,i)=aux(1:3)
             end if
             !
          enddo
       enddo
       !
       end subroutine gaussian
       !
       end module func
       !-----------------------------------------------------------------------
